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ABSTRACT 

The globular cluster population in M87 has decreased measurably through dynam- 
ical evolution caused by relaxation, binary heating and time-dependent tidal pertur- 
bation. For fundamental plane ellipticals in general, cluster populations evolve more 
rapidly in smaller galaxies because of the higher mass density. A simple evolutionary 
model reproduces the observed trend in specific frequency with luminosity for an ini- 
tially constant relationship. 

Fits of theoretically evolved populations to M87 cluster data from McLaughlin et 
al. (1994) show the following: 1) dynamical effects drive evolution in the initial mass 
and space distributions and can account for the large core in the spatial profile as well 
as producing radial-dependence in the mass spectrum; 2) evolution reduces Sn by 50% 
within 16 kpc and 35% within 50 kpc, implying that Sn was initially 26 in this region. 
We estimate that 15% of the 'missing' clusters lie below the detection threshold with 
mass less than 10 5 M Q . 
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1 INTRODUCTION 

Observations of some giant elliptical galaxies reveal globular 
cluster systems which appear more extended than the host 
(Harris 1991). A particularly well-documented example be- 
longs to M87 with a core radius of 7 arc sec and a cluster 
system with a core radius of 1 arc min (McLaughlin 1995). 
However, because a cluster population evolves dynamically 
due to both internal and external processes, the currently 
observed population almost certainly differs from the pri- 
mordial one, complicating the interpretation. 

Researchers have attempted to explain the extended 
core of the M87 cluster distribution as the evolved remnant 
of an initial profile which more closely resembled the light. 
However, neither dynamical friction nor shocking by a com- 
pact nucleus can fully account for this feature. Lauer & Ko- 
rmendy (1986) found that a dynamical friction induced in- 
flow can broaden an initially peaked spatial distribution but 
not at the observed scales. Ostriker, Binney & Saha (1989, 
hereafter OBS) subsequently determined that nuclear tidal 
disruption is viable only if clusters formed exclusively on 
box orbits. 

Another potential mechanism is cluster evaporation 
through dynamical evolution. Recent work in this area 
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demonstrates that evaporative mass loss driven by relax- 
ation and heating due to a time- varying tidal field can lead 
to strong evolution of the Milky Way cluster population in 
a Hubble time (Weinberg 1994, Murali & Weinberg 1996, 
hereafter MW). In this paper, we examine these influences 
on cluster evolution in the dense inner regions of M87 and 
find that they produce the observed flattened profile from a 
peaked initial distribution over a wide range of initial con- 
ditions. Direct estimates of initial conditions using dynam- 
ically evolved parametric models of the spatial distribution 
and cluster mass function indicate that roughly 35% of the 
initial population dissolves or evolves below the detection 
threshold leaving the large core as a result. Furthermore, 
the decay in the size of the cluster population corresponds 
to a decrease in the specific frequency of globular clusters, 
Sn, which denotes the number of clusters per unit galaxian 
luminosity with L measured in units of M v — — 15. 

The high values of Sn found in giant ellipticals have 
become a key point in galaxy formation arguments and sug- 
gest, for example, that the cluster system formed along with 
M87 (e.g. Harris 1991; van den Bergh 1995). Here we show 
that cluster systems decay more rapidly in less luminous 
fundamental plane ellipticals; this leaves larger values of Sn 
in luminous galaxies at the present epoch even if all ellipti- 
cals begin with equal Sn- Our results thus provide at least 
a partial explanation for the observed trend of Sn with L. 
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The plan of the paper is as follows. We summarize our 
choices for the cluster population and the mass model for 
M87 in The assumptions and method for dynamically 
evolving the population is presented in §|j| The main re- 
sults, the statistical comparison of the observed clusters to 
the theoretical models, is presented in This includes an 
exploration of the evolutionary trends, best fit spatial pro- 
files and mass functions, and an inference of the primordial 
population. In §[|, we discuss the discuss the importance of 
the fundamental plane properties on the observed relation 
between specific frequency and luminosity A summary is 
given in §pl 



2 CLUSTER POPULATION 

We assume that the cluster population formed in an ini- 
tial burst approximately HGyr ago. Stellar evolution dom- 
inated cluster evolution for the first Gyr for a Salpeter IMF 
(/? = 2.35) with mi — 0.1 Mq, corresponding to the main 
sequence lifetime of 2 Mq A-star. Our zero-age population 
represents the epoch when, approximately 10 Gyr ago, re- 
laxation, external heating and core collapse heating began 
to drive cluster evolution. 

The fiducial calculations represent zero-age clusters as 
Wo = 5 King models. Comparison calculations using Wo = 7 
clusters show nearly identical evolution over the long time 
scales of interest, in agreement with the results of MW where 
overall evaporation times were found to depend weakly on 
concentration in the range 5 < Wo < 7. We expect similar 
trends in evolution for Wo = 3 clusters (c.f. MW), except in 
high mass, low-eccentricity cases where tidal heating leads 
to rapid disruption. These clusters enhance the destruction 
rate described below, but constitute a very small fraction of 
expected initial populations. 

Each cluster is tidally limited on its orbit in the host. 
While initial cluster densities may differ from the mean den- 
sity required by perigalactic tidal limitation, subsequent evo- 
lution during the first Gyr leads rapidly to tidal truncation 
or disruption. The limiting or tidal radius Rt is uniquely de- 
termined by the cluster mass and orbit. Table [l] summarizes 
the choice of parameters for individual clusters. 

To represent the cluster mass distribution, v(M,r), we 
use pure power laws (e.g. Harris & Pudritz 1994), power 
laws whose exponents have a linear dependence on radius, 
and a Gaussian magnitude distribution (e.g. McLaughlin, 
Harris & Hanes 1994, MHH). Power law mass distributions 
have been proposed on physical grounds by Harris & Pudritz 
(1994) while the Gaussian is commonly used as a convenient 
fitting function for the observed cluster luminosity function. 
To represent the spatial distribution of the cluster popula- 
tion in the primary, we use power law densities with and 
without a core derived from isotropic distribution functions, 
f(E). Orbital isotropy is assumed due to lack of observa- 
tional constraint. 

Adopted models are given by joint distributions 
u(M, r) x f(E) and are summarized in Table |[ The Model 1 
and Model 2 families use power law mass and Gaussian mag- 
nitude distributions respectively. Within each family, succes- 
sive models have additional parameters to explore varying 
core sizes and radial dependence of the mass spectral index. 



Detailed derivation of models from the underlying distribu- 
tion function is given in Appendix 

Finally, we represent the potential of M87 as a singular 
isothermal sphere, with rotation velocity vo = 606 km s - 
(e.g. OBS), velocity dispersion a — 350kms _1 , and assume 
a distance of 16 Mpc (van der Marel 1992). This defines a 
length scale of 77.6 pc per second of arc. Further discussion 
of potential and distance scale is given in Appendix 



3 CLUSTER EVOLUTION 

Competition between internal relaxation and heating due 
to external forcing may dramatically affect a cluster's evo- 
lutionary time scale and survival history. In addition to im- 
pulsive heating of a cluster halo — in a gravitational bulge 
shock, for example — resonances between the cluster's own 
orbital motion and internal stellar trajectories may heat 
cluster stars beyond the limit set by adiabatic invariance 
(Weinberg 1994). For tidally-limited clusters resonant heat- 
ing on low-eccentricity orbits and tidal limitation on high- 
eccentricity orbits drive rapid cluster evolution and evapo- 
ration (see MW for details). The strength of these effects 
motivates this study. 

The evolution of individual clusters includes two-body 
relaxation in the one-dimensional Fokker-Planck approxima- 
tion (e.g. Cohn 1979), external heating due to the time- 
varying tidal field (MW), and a phenomenological binary 
heating term (e.g. Lee et al. 1991). 

We take advantage of the scale-free galaxian profile by 
fixing orbital energy E of all clusters, choosing an initial grid 
of tidally-limited clusters in k = J / J m ax {E) and mass, and 
computing the evolution to complete evaporation. The quan- 
tity J m ax(E) denotes the maximum angular momentum of 
an orbit with energy E. This grid may then be scaled to 
all desired orbital energies. The time evolution of the space 
density for the entire population is then constructed by de- 
termining the phase space distribution at the desired time 
using the evolutionary calculations and projecting appropri- 
ately. 

Although we specifically consider M87, the results apply 
to any elliptical with similar profile. For example, we can 
scale evolution to any fundamental plane elliptical. Because 
the period decreases with mass, the same initial population 
will be more evolved for smaller mass primaries (see §jB| for 
more discussion). 



4 RESULTS 

4.1 Evolution of the core 

An initially peaked cluster distribution develops a flattened 
core through dynamical evolution of individual clusters, as 
shown in Figure [l|. In the dense regions of the inner galaxy, 
rapid mass loss due to relaxation and tidal heating can 
cause complete evaporation of a cluster or drive it below 
the observational limit. Rapid relaxation results from the 
high densities imposed by tidal limitation while tidal heat- 
ing strongly enhances evaporation rates on low-eccentricity 
orbits. The resulting profiles are similar to the profile derived 
in McLaughlin (1995). 

The overall shape of the evolving profile depends on 
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Table 1. Cluster Initial Conditions 



Structural parameters 


Fiducial value 


M total mass 

Wo King concentration parameter 
Rt cluster limiting radius 


10 5 < M c < 5 X 1O 6 M 
Wo = 5, 7 
tidal limitation 


Mass spectral parameters 


/3 mass spectral index: N(m) oc m~ @ 
mi lower mass limit 
m u upper mass limit 


P = 2.35 (Salpeter) 
mi = 0.1 
m u = 2.0 


Orbital parameters 


E orbital energy 

k relative ang. mom.: J/J m ax(E) 


isotropic orbit 
distribution 



Table 2. Population models 



Designation 


p(r) 






v(M) 


Parameters 


Model la 










r), a 


Model lb 








M~(a + Kr) 


r), a, K 


Model lc 


po{r 2 +r 2 )" 


n/2 






r>, r c , a 


Model Id 


po{r 2 c +r 2 Y 


Tj/2 




M - {a + Kr) 


r], r c ,a, K 


Model 2a 
Model 2b 


Por~ v 
po(r 2 +r 2 Y 


>)/2 




-V ) 2 /2al . d y/ dM 
~V a f/2vl . d y/ dM 


ri,V ,a v 
rj, r c , Vo,o- V 




log R (kpc) 



Figure 1. Surface density inferred from observations (open 
squares) compared to Model la for indicated values of mass spec- 
tral index a and initial r -3 (rj = 3) profiles (solid curves). Rapid 
evolutionary rates of low mass clusters produce flatter cores for 
larger a. 

both the initial mass and space distribution of the clusters. 
Consider the following limits. For a fixed spatial profile, a 
distribution rich in low mass clusters evolves rapidly due to 
short evaporation times while for a fixed mass distribution, 
a sharply peaked spatial profile develops a smaller core than 
a shallow profile. Taken together, the two trends produce 



a correlation between the inferred initial mass distribution 
and density profile: a large population of low mass clusters 
can rapidly flatten a steep initial profile while, conversely, a 
large population of high mass clusters allows a flatter initial 
profile to evolve slowly to the same final shape. Observations 
of the cluster mass distribution will distinguish between the 
different initial conditions. 

4.2 Estimates of initial conditions 

The evolved cluster populations described in ^ are com- 
pared to observed cluster dataQ using a maximum likeli- 
hood estimator which combines model and background sur- 
face densities with incompleteness measurements. The back- 
ground surface density is taken to be 6.33 per arcmin 2 with 
a uniform luminosity function (MHH). Point sources lie in 
the region 1.21' < R < 7' (the field edge), centered on M87, 
with apparent limiting magnitude V = 24. We use a mass- 
to-light ratio (M/L)v — 2 to convert luminosity to mass. 
Note that larger M/L will shift the population to higher 
mass, implying less evolution, while smaller M/L will have 
the opposite effect. For R < 1.21', the authors provide 3 
binned points (McLaughlin 1995) . A joint ^ 2 -likelihood es- 
timator is used to include all data points. Details of the 
estimation procedure are provided in Appendix |^. 

We fit both dynamically evolved distributions and un- 
evolved distributions based on the models presented in Ta- 
ble ^. In the case of dynamically evolved models, the quoted 
values represent initial conditions (labeled by 'initial'). In 
the case of unevolved models, the quoted values represent 
the best fit parameters at the present epoch (labeled by 

T The data have been kindly supplied by McLaughlin & Harris 
(1995) 
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Figure 2. 95% and 99% confidences for marginal density in 
Vb and cry for present (dashed) and initial (solid) Model 2a fits. 
Points indicate best-fit values for present and initial models. The 
magnitude limit of the data is indicated as the vertical line. 



'present'). Only models with cores are considered in the 
present epoch fits because coreless models poorly represent 
the data. 

Tables ^ and ^| present the best estimates and their 
variances (cf. Table ^). Comparison of present epoch and 
initial parameter estimates illustrates several expected evo- 
lutionary trends. The core of the distribution grows due to 
the depletion of clusters in the inner regions of the galaxy. 
The power law index a decreases, while the Gaussian mag- 
nitude peak Vb and slope K of the radially-dependent power 
law both increase as a result of the selective evaporation of 
lower mass clusters. The increase in K also indicates that 
depletion occurs primarily in the inner regions. The follow- 
ing cases show specific features of these trends. 

We plot the marginal probability density in Vb and cry 
for Model 2b in Figure ^. The best estimate for Vo decreases 
with time as low mass (and therefore high V) clusters dis- 
appear. However, the assumption of identical initial and 
present epoch values cannot be ruled out since the values 
of Vb are only weakly inconsistent (cf. Table ^) . The lack of 
constraint could result from the shallow magnitude limit of 
the data. Both fits are consistent with distributions which 
peak below the limiting magnitude of V — — 7. 

Model Id suggests that there is radial dependence in 
the mass distribution (Fig. ||). Both present epoch and ini- 
tial fits are inconsistent with a constant mass spectral in- 
dex (K — 0) and indicate that dynamical evolution has en- 
hanced the radial dependence. In the core region, the present 
index a « 1.4. We note that these results conflict with those 
of MHH and McLaughlin & Pudritz (1996), who find no ra- 
dial dependence in the mass distribution. 




1 1.2 1.4 1.6 1.8 2 

a 

Figure 3. 75%, 90%, 95%, and 99% confidences for marginal 
density in a and K for present epoch and initial Model Id fits. 
Points mark best-fit values. The initial spectral index shows mild 
radial dependence increasing from a central value of 1.61 to 1.97 
at 30 kpc in the best-fit case. 



4.3 Comparison of models 

The previous section examined the results of trends in the 
evolution of the cluster population within several model fam- 
ilies. Here, we identify the best overall representation among 
the initial models after 10 Gyr of evolution using a general- 
ized likelihood ratio test (e.g. Martin 1971). 

The final column in Table ^ shows that the most general 
model gives a better estimate in each case, as expected. In 
particular, Model lc can be rejected in favor Model Id, a 
result which is consistent with the confidence surfaces for 
the radially-dependent mass spectrum plotted in Figure ^. 

Model 2b generalizes Model 2a by introducing arbitrary 
initial core size. A finite core does provide a better estimate 
but zero core (Model 2a) cannot be rejected. Figure ^ com- 
pares the surface density profiles of these two models. The 
Model 2b fit falls below the observed profile at small radii 
due to the shallow initial core. However, the binned surface 
density points have relatively low weight in the full data set. 
The good fit of Model 2a to the inner data points suggests 
that a more peaked initial distribution may provide the op- 
timal fit. The use of individual cluster counts in this region 
should help provide the necessary constraint. The deviation 
of the present epoch fit from the data further suggests that 
evolution plays an important role in shaping the profile. 

Finally, we compare the most general power law mass 
function model, Model Id, with the most general Gaussian 
magnitude model, Model 2b, by constructing a linear combi- 
nation of both spaces and searching for the global maximum. 
The maximum occurs at the best-fit parameters for Model 
2b: the Gaussian magnitude distribution describes the data 
significantly better than any power law mass distribution. As 
discussed above, the Gaussian may be poorly constrained by 
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Table 3. Model 1 fits 



Table 4. Model 2 fits 



Epoch 


V 


O" 17 


r c 






K 


CK 


-log L 


Model la 


initial 


2.76 


0.04 




1.95 


0.03 






69502.9 


Model lb 


initial 


2.66 


0.06 




1.68 


0.09 


0.01 


0.004 


69501.3 


Model lc 


present 


3.09 


0.13 


7.30 


0.79 1.73 


0.04 






69500.6 


initial 


3.03 


0.12 


5.67 


0.94 1.93 


0.04 






69499.8 


Model Id 


present 


3.06 


0.13 


7.34 


0.81 1.30 


0.08 


0.028 


0.002 


69480.9 


initial 


3.13 


0.10 


5.70 


0.40 1.61 


0.10 


0.014 


0.004 


69492.8 


Epoch 


'/ 


O it 


r c 


o> c Vq 


°v 


ay 


0V V 


-log L 


Model 2a 


initial 


2.77 


0.05 




-7.11 


0.17 


1.16 


0.08 


69487.2 


Model 2b 


present 


3.10 


0.13 


7.32 


0.80 -7.33 


0.10 


1.08 


0.07 


69468.7 



initial 3.14 0.12 5.14 0.77 -7.07 0.17 1.19 0.08 69485.3 
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Figure 4. Comparison of surface density profiles of Models 2a 
and 2b (solid) with binned data. Initial profiles (dotted) and 
present epoch Model 2b fits (dashed) are also plotted. Both the 
model with initial core and the present epoch fit deviate from the 
data in the inner region. 



the limiting magnitude of the data. However, examination 
of the estimated functions shows that the power law is more 
peaked at low mass than the Gaussian both initially and fi- 
nally, while the Gaussian has more weight at high mass. Thes 
differences in shape also lead to the statistical preference of 
the Gaussian luminosity function over the spatially constant 
power law and its radially-dependent generalization. 



Table 5. Evolved model comparisons 
test — 21nA 1 ' lA accept reject confidence 



lc-ld 
2a-2b 
ld-2b 



14.0 
3.8 
15.0 



V 



V 
V 



99% 
60% 
96% 



likcliho 



Table ^ summarizes the conclusions of the tests com- 
paring the initial conditions and lists the likelihood ratios 
and confidence values. In the first two cases, the likelihood 
ratios follow directly from Tables ^| and ^. In summary, we 
find that the initial Gaussian magnitude distribution best 
describes the data, but that we cannot distinguish between 
singular density profiles and densities with core. Both con- 
clusions may result from insufficient data. 



4.4 Evolution of the initial population 

From the derived initial conditions, we plot the projected 
cumulative distribution of clusters initially within 50 kpc 
or 6.5R e (Figure |E|). The initial population in this region 
is about 7250 for Model 2a, a factor of 1.6 larger than the 
presently observed population of 4500. For Model 2b, the 
initial population size is smaller, but still in excess of 6000. 
Evolved clusters will also be found at masses below the ob- 
servational limit. Within 16 kpc, about 500 evolved clusters 
are expected in the range 24 > V > 24.5 (10 5 > M > 
7 x 10 4 Af©). All other clusters in this region initially in this 
mass range will have evaporated. 

This implies that the specific frequency, Sjv, evolves 
with time. Using the ratio of final to initial surface densities 
from the models, we derive the run in initial specific fre- 
quency at radius R, Sn(R) (Figure^). As expected, deple- 
tion in the inner regions dominates; however, ~ 10% change 
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Figure 5. Estimated cumulative distribution of clusters for pro- 
jected radius R < 50 kpc (6.5R e ). Initial distributions for Model 
2a (dashed) and Model 2b (solid) are shown for V < 24 (upper 
pair). Final distributions are shown for V < 24 (middle pair) and 
for 24 > V > 24.5 (bottom pair). Clusters with 24 > V > 24.5 
began with V < 24. Approximately 50% of the initial clusters 
vanish within (16 kpc) 2R e and 35% within 6.5-Re for Model 2a. 



Figure 6. The ratio of final to initial surface density (top) in 
Models 2a (dashed) and 2b (solid) and the run in initial specific 
frequency at R (Sff(R)) for each model (bottom) derived from 
the observed values given by MHH (dotted). Evolution reduces 
S N by 35% in Model 2a. 



in Sn (R) occurs even out to 50 kpc due to the rapid evolu- 
tion of low-mass clusters. 

Model 2a in Figure |^ shows a 35% change in total Sn 
within 50 kpc due to depletion. The observed total value 
Sat ~ 17 in this region (MHH) implies an initial value of 
26.5. Thus even the enormous observed value of Sn has 
diminished significantly due to evolution (this neglects the 
intrinsic evolution in galaxy luminosity due to stellar evolu- 
tion). The time evolution of the total Sn is shown directly 
in Figure The decay of the cluster population is approx- 
imately exponential in time with e-folding times of 20 Gyr 
and 40 Gyr for measurements within 16 and 50 kpc, respec- 
tively. 

Our comparison applies to clusters in specified mass and 
radial ranges in a galaxy. Quoted specific frequency values 
are extrapolated over unobserved ranges from data taken 
within such limits. In the case of M87, MHH find S N = 17.7 
directly from the observations, while extrapolation yields a 
total Sn = 14.4 for the whole galaxy. Their extrapolation 
of the luminosity function over all magnitudes yields a total 
correction of 2.2 to cluster counts in the observed magnitude 
range. For R < 42 kpc, using this correction they estimate 
a total ~ 9400 clusters (3729 observed in galaxy+cD enve- 
lope, 500 in core, x2.2) which is 70% of the total number of 
clusters estimated over all radii. The estimated initial dis- 
tribution has a larger corr ectio n because it is weighted more 
heavily to low mass (c.f. §4.2). Since most of these clusters 
completely evaporate, this implies even greater evolution in 
specific frequency than derived here. 
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Figure 7. The fraction of clusters remaining within the indicated 
projected radius as a function of time. Abscissae indicate time 
units which correspond to fundamental plane scaling for galaxies 
with indicated masses. Top axis gives M87 scaling (other axes will 
be discussed in the next section). Initial Sn is ~ 21 measured 
within 16 kpc and 26.5 measured within 50 kpc and decays with 
the cluster population. 
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5 DISCUSSION 



Our conclusions ignore the possibility that recent merger 
and accretion events have have strongly contaminated the 
initial cluster population. Merging of gas-rich galaxies is ex- 
pected to produce clusters with strong central concentration 
(Zepf & Ashman 1993; Mihos & Hernquist 1994), but the 
large core of the distribution itself argues against any recent 
merger which has produced significant numbers of young 
clusters. Thus either cluster-producing mergers in M87 have 
occurred long in the past or not at all. 

The recent addition of clusters through satellite accre- 
tion is also unlikely to account for a significant fraction of 
the observed population. For example, accretion of Milky 
Way-type spirals can only account for about 25% of the ob- 
served population, assuming that the 4L, luminosity of M87 
(R < 40kpc) comes entirely from satellites (Lauer 1988). 
Assuming that material is stripped when the mean density 
of the primary exceeds that in the satellite, accretion will 
deposit clusters in regions of mean density similar to that in 
the original environment of the accreted galaxy. The clus- 
ters, then, will remain roughly tidally truncated. Their new 
orbits depend on the orbit of the dissolving satellite, but 
otherwise evolution should be similar and the accreted pop- 
ulation may appear coeval regardless of the time of accre- 
tion. 

The high specific frequency of globular clusters in M87 
appears to indicate the exceptional conditions governing the 
formation and evolution of cD galaxies relative to other ellip- 
ticals. The observation that specific frequency increases with 
galaxy luminosity suggests that galaxy formation was not an 
intrinsically hierarchical and homologous process (e.g. San- 
tiago & Djorgovski 1993). 

However, the current results indicate that density dif- 
ferences between galaxies will lead to differences in the evo- 
lution of intrinsic cluster populations. Since M87 lies ap- 
proximately on the fundamental plane for elliptical galaxies, 
we can investigate differences in environment-driven cluster 
evolution by scaling our results to other ellipticals. We as- 
sume that the initial profile of the cluster population de- 
rived above, when scaled homologously, describes the initial 
population in any elliptical. The dynamical time scale for a 
tidally truncated cluster with constant mass and spatial pro- 
file is then determined by the mean density of the galaxy: 
r ~ M- 1/2 R 3/2 . This yields r ~ M 13 , where = 0.4 for 
the fundamental plane scaling given in Faber et al. (1987), 
= 0.26 for the Djorgovski & Davis (1987) results, and 
— 0.6 for a more recent set of parameters from Faber 
(1995; see also Pahre et al. 1995). This relation implies that 
clusters evolve and are depleted more rapidly in smaller, 
low luminosity ellipticals than in massive, high luminosity 
ellipticals. 

To demonstrate the importance of this effect, we com- 
bine the approximat e ex ponential decay rate of the cluster 
population found in §4.4 with this scaling relation. This gives 



an expression for the number of clusters remaining in an ini- 
tially coeval population belonging to a galaxy of luminosity 
L at time T: 



N c i(L,T) = jVo(L)e- T/T ° (AW/M) " 
= N e~ T/To{LM87/L)1 ' 2 ' L ' 3 , 



@-0A (solid) 
log N,-4.00±.05 
7=1.27±.13 




jS=0.6 (dotted) 
log N.=3.97±,05 
_ 7=1. 03 = . 13 



-24 
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-22 
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-20 
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Figure 8. Evolved population size as a function of galaxy lumi- 
nosity compared to data from Harris (1996) for fundamental plane 
parameter /3 = 0.4 (solid) and /3 = 0.65 (dotted) . The present- 
day, unevolvcd fit (dashed) has logiV, = 3.89±.05,7 = 1.46±.13. 
For /3 = 0.65, the initial value of 7 is consistent with specific fre- 
quency which is independent of luminosity 



where to ~ 40Gyr for M87, No(L) is the initial distribution 
of cluster population sizes as a function of galaxy luminosity 
and M/L ~ L 1,24 throughout. 

Using a power-law distribution No(L) — N,(L/ Lmsi)' 1 , 
we compare the model curve N c i(L,T) to an observed sam- 
ple of cluster populations in galaxies (Harris 1996). The re- 
sulting models show qualitative agreement with the data 
(Figure ^ , falling off at low luminosity more rapidly with in- 
creasing due to the more rapid rate of evolution. From this 
we conclude that the observed discrepancies in specific fre- 
quency which correlate with galaxy luminosity were smaller 
in the past. This is reflected in the smaller exponent 7 in 
N c i oc L~' predicted initially compared to the present-day 
estimate. 

To conclude, we suggest that cluster evolution may ac- 
count at least in part for the specific frequency problem. 
Because clusters evolve more rapidly in dense environments 
and because ellipticals are typically denser at low luminosity 
that at high luminosity, cluster populations diminish more 
rapidly in low luminosity galaxies. Specific frequencies will 
therefore correlate more strongly with luminosity at later 
times. 



6 SUMMARY 

We have investigated the evolution of the M87 globular clus- 
ter system. The results of this study have enabled us to ex- 
amine the broader question of population evolution in fun- 
damental plane elliptical galaxies. Our main conclusions fol- 
low: 

(i) The loss of globular clusters through relaxation and 
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tidally-induced evaporation accounts for the large core and 
shallow profile in cluster number distribution compared to 
the light distribution in M87. 

(ii) Evolution produces a radial dependence in the 
present-day mass spectrum of clusters such that higher mass 
clusters predominate in the inner regions. The models also 
indicate an initial radial dependence in the mass spectrum. 

(iii) Likelihood ratio tests reject an initial power law in 
favor of an initial Gaussian but do not rule out the possibility 
of an initial core in the profile. 

(iv) The best-fit model for M87 has an initial population 
of 7.25 x 10 3 clusters with projected radius ii < 50 kpc, 
about 60% more than is currently observed. Roughly 14% 
of the initial population are now objects of slightly less than 
10 5 M©; dynamical evolution can strongly modify the spe- 
cific frequency of globular clusters. 

(v) Scaling the calculations to fundamental plane ellip- 
tical galaxies indicates that cluster evolution in the differ- 
ing environments qualitatively accounts for the trend in ob- 
served population number versus galaxy luminosity. Smaller 
galaxies tend to have high densities and thus more rapid 
evolutionary time scales, so their cluster populations tend 
to diminish more rapidly. 

Some of these inferences, especially (iii), may be biased 
by the lack of detailed data in the inner galaxy. Point source 
data for r < 1' will lead to stronger constraints on the size 
of the initial core and a deeper survey will pick up low lumi- 
nosity objects and further constrain the luminosity function. 
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APPENDIX A: CLUSTER DISTRIBUTION 
FUNCTIONS 

The initial cluster population is represented by joint dis- 
tributions of the phase space density, f(E) and the mass 
spectrum v(M,r): 



i>(M,r,E,J) 



dN 



x f(E,r)u(M,r), 



(Al) 



dMdEdJ 

where ip(M, r, E, J) is the number of stars per unit mass per 
unit energy per unit angular momentum at a fixed point in 
space. 

For the initial orbit distribution in coreless models, we 
generalize the isothermal distribution employed by OBS 

Po 



HE) 



-E/o 



(27^)3/2 



(A2) 



with o"o = t]~ 1 Vq. Note that rj^ 1 — | for the isothermal 
sphere and = § for the fiducial model employed by OBS. 
The background gravitational potential of M87 is taken to 
be a singular isothermal sphere, = «o mr i indepen- 

dent of the cluster distr ibut ion. The space number density 
of clusters for equation (A2) is then 



n(r) = p r n . 

We generalize this profile to include a core: 
n(r) = po{r 2 + r 2 y v/2 



(A3) 



(A4) 



where r c is the core radius of the system. The isotropic clus- 
ter distribution function follows from integral inversion (e.g. 
Binney & Tremaine 1987). 
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Figure Bl. Comparison of singular isothermal sphere (SIS), soft- 
ened isothermal sphere used by OBS, and the luminosity profile 
determined by van der Marel (1994). The isothermal models agree 
with the observed profile at small radii and continue to rise lin- 
early with the dark matter halo. 



We consider three initial distributions of cluster masses: 
1) a Gaussian distribution of initial magnitudes, V , which is 
everywhere constant in space (e.g MHH); 2) a power law dis- 
tribution of mass, M, which is everywhere constant in space 
(e.g. Harris & Pudritz 1994) ; and 3) a power law distribution 
of mass with a radially-dependent spectral index. 

The Gaussian distribution of initial magnitudes V de- 
fines the mass spectrum 



v(M) oc e 



(- 



-1-0)2 



dV 
dM' 



(A5) 



where the transformation to mass is effected by the Jaco- 
bian, dV/dM. This distribution is characterized by two pa- 
rameters: Vo and ov, the mean and dispersion of magni- 
tudes. 

The simple power law mass distribution 



v(M) oc M~ 



(A6) 



depends only on the mass spectral index, a. We define the 
following radially dependent distribution 

-(a + Kr) 



v{M,r) oc M~ 



(A7) 



whose spectral index has the central value a and varies lin- 
early with radius. 



second of arc. The more recent estimate by Elson & Santiago 
(1996) based on the Cepheid calibration of Freedman et al. 
(1994) gives a length scale of 87 pc per second of arc which is 
closer to the latter estimate but not strongly discrepant from 
the first. OBS apparently adopted the length scale of 132 pc 
per second of arc, corresponding to Ho = SOkms" 1 Mpc -1 . 

For the singular isothermal sphere (SIS), <!>(r) = Vq lnr, 
we adopt vo = 606 km s" 1 (OBS). This defines a velocity 
dispersion of a = 350 km s -1 for an initial distribution of 



luminous matter with r/ 1 = i. Figure Bl compares the 



mass of the SIS with the softened isothermal sphere used by 
OBS, and van der Marel's (1994) estimate derived using the 
luminosity profile and dynamical modeling. 

The results presented here do not significantly depend 
on choice of singular or softened potential. The 10% mass 
difference at 3 kpc leads to 5% faster evolution in the SIS 



since the dynamical time scale goes as p 



-1/2 



in tidally- 



limited clusters. This is well within the 4.7 kpc core radius of 
the cluster system determined by McLaughlin (1995). Con- 
versely, the tidal field in the the slowly rising core region of 
the softened model is closer to that in a Keplerian potential 
which is stronger and so balances the larger mass of the SIS. 



APPENDIX C: MAXIMUM LIKELIHOOD 
ESTIMATION OF MODEL PARAMETERS 

A joint x 2 - max i m um likelihood estimator is used in §^ 
to fit the dynamical models to V-band photometric data 
and binned surface density data of the M87 cluster system 
(MHH). THe expected surface density profiles are the sum 
of the dynamically evolved model surface density, S(r,v;9), 
derived from distributions in and background surface 
density, ao(V), multiplied by the incompleteness factor, 
f(x, y, V) which represents the probability of detecting a 
cluster of given magnitude at a particular location in the 
field. The likelihood statistic 



L = [| /(%• , yj , Vj ) [S( rj , Vt ; 9) + a ( V 3 )] , 



(CI) 



defines the posterior probability of the data given the model. 
The x 2 statistic 



nl i ( s^s{ r ,e) Y 



(C2) 



defines the posterior probability of the binned data given 
the model. The total posterior probability is then 



P = P v 2 x L 



(C3) 



and the parameters which maximize the posterior probabil- 
ity of the data are the best estimates. 



APPENDIX B: GENERALIZED ISOTHERMAL 
SPHERE 

We adopt a distance of 16 Mpc to M87 which corresponds to 
Ho = 81.5 kins -1 Mpc -1 (van der Marel 1994). This defines 
a length scale of 77.6 pc per second of arc. Lauer & Kor- 
mendy (1986) adopted Ho = 75kms _1 Mpc" 1 , giving a dis- 
tance of 17.4 Mpc and defining a length scale of 84.1 pc per 



